#########################################################################
#########################################################################
########## Political Considerations in Nonpolitical Decisions: ##########
##########      A Conjoint Analysis of Roommate Choice         ##########
########## 		    	   Richard Shafranek                       ##########
########## 		 	 	     3.February.2019                         ##########		
#########################################################################
#########################################################################
library(tidyverse)
library(readr)
library(ggplot2)
library(dplyr)
library(cjoint)
library(multiwayvcov)
library(miceadds)
library(MASS)

#########################################################################
# read the data

conjoint <- read.csv("conjoint.csv")

#########################################################################
#relevel the factors to use the reference levels we want
conjoint <- within(conjoint, traits_bedtime <- relevel(traits_bedtime, ref = "11 pm"))
conjoint <- within(conjoint, traits_clean <- relevel(traits_clean, ref = "Somewhat clean and tidy"))
conjoint <- within(conjoint, traits_race <- relevel(traits_race, ref = "White"))
conjoint <- within(conjoint, traits_lgbt <- relevel(traits_lgbt, ref = "Straight"))
conjoint <- within(conjoint, traits_religion <- relevel(traits_religion, ref = "Nonreligious"))
conjoint <- within(conjoint, traits_pid <- relevel(traits_pid, ref = "Independent"))
#create a non-releveled version of political interest for use in interaction models
conjoint$traits_polinterest_n <- conjoint$traits_polinterest
conjoint <- within(conjoint, traits_polinterest <- relevel(traits_polinterest, ref = "Somewhat interested in politics"))
conjoint <- within(conjoint, traits_hobby <- relevel(traits_hobby, ref = "Shopping"))
conjoint <- within(conjoint, traits_value <- relevel(traits_value, ref = "Trying new things"))

#########################################################################
# AMCEs
# we need to do some subsetting, because A) we're only interested in partisans here (independents don't have an "outparty"), and B) AMCEs only work with observations without missing data, so we need to drop incomplete cases
conjoint_complete <- drop_na(conjoint, c(rating, choice, traits_pid_relative))
conjoint_complete_partisan <- conjoint_complete[conjoint_complete$pid3_withlean!="Independent"&conjoint_complete$pid3_withlean!="Other party (please specify)",]

#we need to use makeDesign since the attribute distribution is not uniform (due to the weighted levels of LGBT)
#first, list the attributes
attribute_list <- list()
attribute_list[["traits_race"]] <- c("White", "Asian American", "Black", "Hispanic")
attribute_list[["traits_lgbt"]] <- c("Straight", "LGBT")
attribute_list[["traits_religion"]] <- c("Catholic", "Evangelical Christian", "Nonreligious", "Jewish")
attribute_list[["traits_pid_relative"]] <-c("Inparty","Independent", "Outparty")
#attribute_list[["traits_pid"]] <-c("Democrat","Independent", "Republican")
attribute_list[["traits_polinterest"]] <-c("Not at all interested in politics", "Somewhat interested in politics", "Very interested in politics")
attribute_list[["traits_clean"]] <- c("Not at all clean and tidy", "Somewhat clean and tidy", "Very clean and tidy")
attribute_list[["traits_music"]] <-c("Classic rock", "Country music", "Electronic", "Hip-hop", "Jazz")
attribute_list[["traits_hobby"]] <-c("Shopping", "Cars & auto mechanics", "Doing yoga", "Hunting and fishing",
                                     "Playing golf", "Swimming", "Theatre/performing arts", "Visiting farmers markets",
                                     "Watching foreign films", "Watching sports")
attribute_list[["traits_social"]] <- c("Likes to go to parties on weekends", "Likes to stay in on weekends")
attribute_list[["traits_value"]] <- c("Trying new things", "Following rules and behaving properly", 
                                      "Helping others around them", "Respecting traditions", "Treating others fairly")
attribute_list[["traits_bedtime"]] <-c("11 pm", "1 am", "9 pm")

marginal_weights <- list()
marginal_weights[["traits_race"]] <- rep(1/length(attribute_list[["traits_race"]]), length(attribute_list[["traits_race"]]))
marginal_weights[["traits_lgbt"]] <- c(4/5, 1/5)
marginal_weights[["traits_religion"]] <- rep(1/length(attribute_list[["traits_religion"]]), length(attribute_list[["traits_religion"]]))
marginal_weights[["traits_pid_relative"]] <- rep(1/length(attribute_list[["traits_pid_relative"]]), length(attribute_list[["traits_pid_relative"]]))
marginal_weights[["traits_polinterest"]] <- rep(1/length(attribute_list[["traits_polinterest"]]), length(attribute_list[["traits_polinterest"]]))
marginal_weights[["traits_clean"]] <- rep(1/length(attribute_list[["traits_clean"]]), length(attribute_list[["traits_clean"]]))
marginal_weights[["traits_music"]] <- rep(1/length(attribute_list[["traits_music"]]), length(attribute_list[["traits_music"]]))
marginal_weights[["traits_hobby"]] <- rep(1/length(attribute_list[["traits_hobby"]]), length(attribute_list[["traits_hobby"]]))
marginal_weights[["traits_social"]] <- rep(1/length(attribute_list[["traits_social"]]), length(attribute_list[["traits_social"]]))
marginal_weights[["traits_value"]] <- rep(1/length(attribute_list[["traits_value"]]), length(attribute_list[["traits_value"]]))
marginal_weights[["traits_bedtime"]] <- rep(1/length(attribute_list[["traits_bedtime"]]), length(attribute_list[["traits_bedtime"]]))

roommate_design <- makeDesign(type="constraints", J=NULL, filename=NULL, attribute.levels=attribute_list, constraints=NULL,
                              level.probs=marginal_weights, tol=1e-14)

#now, run the analyses
#this becomes Figure 2
fig2_results <- amce(rating ~ traits_race + traits_lgbt + traits_religion + traits_pid_relative + traits_polinterest
                     + traits_clean + traits_music + traits_hobby + traits_social + traits_value + traits_bedtime, data = conjoint_complete_partisan, design = roommate_design, respondent.varying = NULL, subset = NULL,
                     respondent.id = "ResponseId", cluster = TRUE, na.ignore = TRUE, weights = NULL, baselines = NULL)

#plot:
plot(fig2_results, main = "", xlab = "Change in E[Y] (with 90 percent confidence intervals)", ci = 0.90,
     colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7",
                "#994F00", "#000000", "#5D3A9B"), xlim = NULL, breaks = NULL,
     labels = NULL, attribute.names = c("Bedtime", "Cleanliness","Hobby", "Sexual orientation",
                                        "Music", "Partisanship", "Political interest", "Race",
                                        "Religion", "Social preferences", "Value"), level.names = NULL,
     label.baseline = TRUE, text.size = 10, text.color = "black",
     point.size = 0.5, dodge.size = 0.9, plot.theme = NULL,
     plot.display = "all",
     facet.names = NULL, facet.levels = NULL,
     group.order = NULL,font.family = NULL) 

#########################################################################
#match model with OLS and clustered standard errors 
#this becomes Figure 3

#full sample (use for Figure 3, Table 5, and Table A3):
match_cluster_full <- lm.cluster(conjoint, rating ~ traits_pid_match + traits_interest_match + traits_social_match + traits_lgbt_match + traits_religion_match + 
                                   traits_bedtime_match + traits_value_match + traits_race_match + traits_clean_match
                                 + traits_hobbies_match + traits_music_match, cluster=conjoint$ResponseId)

#create Figure 3
t3 <- cbind(coef = summary(match_cluster_full)[2:12], conf_minus = confint(match_cluster_full, level=0.90)[2:12,1], conf_plus = confint(match_cluster_full, level=0.90)[2:12,2], order=c(1:11)) %>% data.frame() 

t3$names <- c("Partisan affiliation match", "Political interest match", "Social match", "Sexual orientation match", "Religion match", 
           "Bedtime match", "Value match", "Race match", "Cleanliness match", "Hobby match", "Music match")
t3$names <- factor(t3$names, levels = t3$names[order(t3$coef, decreasing=FALSE)])

fig3 <- ggplot(t3, aes()) + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + geom_linerange(aes(x = t3$names, ymin = t3$conf_minus, ymax = t3$conf_plus, colour="red"),lwd = 1, position = position_dodge(width = 1)) + geom_point(aes(x = t3$names, y=t3$coef, colour="red")) + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + coord_flip() + theme_bw() + theme(axis.title.y = element_blank(), legend.position = "none", axis.text.y=element_text(hjust=0)) + ylab("Coefficient (with 90 percent confidence intervals)")
fig3	


#########################################################################
#create Figure 4
#subset and relevel
conjoint <- within(conjoint, traits_polinterest_plot <- relevel(traits_polinterest, ref = "Not at all interested in politics"))

#here we only want subjects who are sorted into inparty/outparty on this variable (i.e., exclude non-partisans)
conjoint$traits_pid_relative_only <- ifelse(conjoint$traits_pid_relative=="Inparty", "Inparty", 
                                                  ifelse(conjoint$traits_pid_relative=="Outparty", "Outparty", NA))

#we also only want to look at the "not at all" and "very" levels of political interest (for the sake of this plot)
conjoint$traits_polinterest_plot <- ifelse(conjoint$traits_polinterest_plot=="Not at all interested in politics", "Not at all interested in politics",
                                                 ifelse(conjoint$traits_polinterest_plot=="Very interested in politics", "Very interested in politics", NA))
conjoint_clean_plot <- subset(conjoint, traits_pid_relative_only=="Inparty" | traits_pid_relative_only=="Outparty")
conjoint_clean_plot <- subset(conjoint_clean_plot, traits_polinterest_plot=="Not at all interested in politics" | traits_polinterest_plot=="Very interested in politics")

conjoint_clean_plot %>% 
  group_by(traits_pid_relative_only, traits_polinterest_plot) %>% 
  summarise(rating_groups = mean(rating, na.rm=TRUE)) -> interest_plot

#get confidence intervals with t.test
interest_plot$ymin <- c(t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Inparty"&conjoint_clean_plot$traits_polinterest_plot=="Not at all interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[1],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Inparty"&conjoint_clean_plot$traits_polinterest_plot=="Very interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[1],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Outparty"&conjoint_clean_plot$traits_polinterest_plot=="Not at all interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[1],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Outparty"&conjoint_clean_plot$traits_polinterest_plot=="Very interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[1])
interest_plot$ymax <- c(t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Inparty"&conjoint_clean_plot$traits_polinterest_plot=="Not at all interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[2],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Inparty"&conjoint_clean_plot$traits_polinterest_plot=="Very interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[2],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Outparty"&conjoint_clean_plot$traits_polinterest_plot=="Not at all interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[2],
                        t.test(conjoint_clean_plot$rating[conjoint_clean_plot$traits_pid_relative_only=="Outparty"&conjoint_clean_plot$traits_polinterest_plot=="Very interested in politics"], na.rm=TRUE, conf.level=.90)$conf.int[2])

#plot
interest_plot %>% 
  ggplot() +
  aes(x = traits_polinterest_plot, y = rating_groups, color = traits_pid_relative_only, group=traits_pid_relative_only) +
  geom_line(aes(group = traits_pid_relative_only, linetype=traits_pid_relative_only), size=1) +
  scale_linetype_manual(values=c("solid", "dotdash")) +
  geom_point(aes(shape=traits_pid_relative_only), size = 4) + 
  geom_errorbar(aes(ymax=ymax, ymin=ymin), width=.025) -> interaction_plot

interaction_plot$labels[3] <- "Roommate \npartisanship"
interaction_plot$labels[5] <- "Roommate \npartisanship"
interaction_plot$labels[6] <- "Roommate \npartisanship"
interaction_plot$labels[2] <- "Rating"
interaction_plot$labels[1] <- "Roommate's level of political interest"
interaction_plot <- interaction_plot + theme_bw() + theme(text=element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
interaction_plot

#########################################################################
#predict preferences given particular values (Table 4)

#generate a stereotypic republican profile
stereotypic_republican <- data.frame(traits_race="White", traits_lgbt="Straight", traits_religion="Evangelical Christian", 		
                                     traits_pid="Republican", 
                                     traits_polinterest="Somewhat interested in politics", traits_clean="Somewhat clean and tidy",
                                     traits_music="Country music", traits_hobby="Hunting and fishing", traits_social="Likes to go to parties on weekends",
                                     traits_value="Respecting traditions", traits_bedtime="11 pm")

#generate a counter-stereotypic republican profile
counter_republican <- data.frame(traits_race="Black", traits_lgbt="LGBT", traits_religion="Nonreligious", traits_pid="Republican", 
                                 traits_polinterest="Somewhat interested in politics", traits_clean="Somewhat clean and tidy",
                                 traits_music="Hip-hop", traits_hobby="Theatre/performing arts", traits_social="Likes to go to parties on weekends",
                                 traits_value="Treating others fairly", traits_bedtime="11 pm")

#generate a stereotypic democratic profile
stereotypic_dem <- data.frame(traits_race="Black", traits_lgbt="LGBT", traits_religion="Nonreligious", traits_pid="Democrat", 
                              traits_polinterest="Somewhat interested in politics", traits_clean="Somewhat clean and tidy",
                              traits_music="Hip-hop", traits_hobby="Theatre/performing arts", traits_social="Likes to go to parties on weekends",
                              traits_value="Treating others fairly", traits_bedtime="11 pm")

#generate a counter-stereotypic democratic profile
counter_dem <- data.frame(traits_race="White", traits_lgbt="Straight", traits_religion="Evangelical Christian", traits_pid="Democrat", 
                          traits_polinterest="Somewhat interested in politics", traits_clean="Somewhat clean and tidy",
                          traits_music="Country music", traits_hobby="Hunting and fishing", traits_social="Likes to go to parties on weekends",
                          traits_value="Respecting traditions", traits_bedtime="11 pm")

#subset the data and run some models based on respondent partisanship
dem_only <- subset(conjoint, pid3_withlean=="Democrat")
rep_only <- subset(conjoint, pid3_withlean=="Republican")

dem_model <- lm(rating ~ traits_race + traits_lgbt + traits_religion + traits_pid * traits_polinterest
                + traits_clean + traits_music + traits_hobby + traits_social + traits_value + traits_bedtime, data=dem_only)
rep_model <- lm(rating ~ traits_race + traits_lgbt + traits_religion + traits_pid * traits_polinterest
                + traits_clean + traits_music + traits_hobby + traits_social + traits_value + traits_bedtime, data=rep_only)

#get predicted ratings by party for each profile
predict(dem_model, stereotypic_republican)
predict(dem_model, counter_republican)
predict(dem_model, stereotypic_dem)
predict(dem_model, counter_dem)
predict(rep_model, stereotypic_republican)
predict(rep_model, counter_republican)
predict(rep_model, stereotypic_dem)
predict(rep_model, counter_dem)

#########################################################################
#trait importance (self-report) ratings for Table 5
summary(conjoint$cleanliness_rank)
summary(conjoint$bedtime_rank)
summary(conjoint$values_rank)
summary(conjoint$social_rank)
summary(conjoint$hobbies_rank)
summary(conjoint$political_rank)
summary(conjoint$music_rank)
summary(conjoint$religious_rank)
summary(conjoint$polinterest_rank)
summary(conjoint$lgbt_rank)
summary(conjoint$race_rank)

#########################################################################
#Table A1
table(conjoint$year_in_school)
table(conjoint$race)
table(conjoint$pid3)
summary(conjoint$age)
table(conjoint$gender)
summary(conjoint$r_lgbt)

#########################################################################
#Table A2
summary(fig2_results)

#########################################################################
#Table A3
summary(match_cluster_full)

#########################################################################
#Table A4
tablea4 <- amce(choice ~ traits_race + traits_lgbt + traits_religion + traits_pid_relative + traits_polinterest
                            + traits_clean + traits_music + traits_hobby + traits_social + traits_value + traits_bedtime, data = conjoint_complete_partisan, design = roommate_design, respondent.varying = NULL, subset = NULL,
                            respondent.id = "ResponseId", cluster = TRUE, na.ignore = TRUE, weights = NULL, baselines = NULL)
summary(tablea4)

#########################################################################
#Table A5
tablea5 <- glm(choice ~ traits_pid_match + traits_interest_match + traits_social_match + traits_lgbt_match + traits_religion_match + 
                     traits_bedtime_match + traits_value_match + traits_race_match + traits_clean_match
                   + traits_hobbies_match + traits_music_match, data=conjoint, family="binomial")
summary(tablea5)

